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ABSTRACT 


In  the  present  study,  a  three-dimensional  differential  field  model  based 
on  general  orthogonal  coordinate  systems  is  developed.  The  model,  which  is 
essentially  designed  for  closed  spaces,  also  includes  the  physical  effects  of 
turbulence,  strong  buoyancy,  full  compressibility,  pressure  rise  due  to  fire 
loading,  surface-surface  and  surface-flame  radiation  exchange,  and  heat  losses 
through  the  wall.  It  is  based  on  a  control-volume  staggered-cel 1  finite- 
difference  approach  with  primitive  variables.  Results  of  numerical  calcula¬ 
tions  based  on  the  field  model  are  compared  with  test  data  for  a  methanol  fire 
in  the  NRL  FIRE  I  test  facility  which  is  in  the  form  of  a  closed  pressure 
vessel.  Reasonable  comparisons  of  the  resulting  pressure  and  temperatures  at 
several  locations  have  been  obtained.  Also  shown  are  the  detailed  velocity 
and  temperature  fields  inside  the  vessel  at  different  time  instants  after  the 
commencement  of  the  fire. 
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INTRODUCTION 

It  is  now  generally  recognized  that  fire  hazards  in  rooms,  compartments, 
corridors,  passageways,  and  other  confined  spaces  can  only  be  minimized  and 
mitigated  by  a  multi-prong  approach  involving  different  strategies.  Among  them 
is  fire  modeling  which  deals  with  the  development  of  mathematical  models  which 
describe  the  physical  and  chem-ical  processes  in  the  spread  of  fire  in  spaces  as 
a  function  of  the  ignition  source,  space  geometry,  and  material  content.  Once 
validated  by  experiments  in  small-scale  laboratory  tests  and/or  full-scale  fire 
tests,  these  models  can  be  used  as  computer-based  simulation  models  to  determine 
the  effects  of  the  significant  parameters  on  the  fire-spread  phenomena.  The 
results  can  then  be  utilized  for  developing  fire  mitigation  measures  as  well  as 
for  providing  a  rational  basis  for  post-fire  investigations.  One  chief  advan¬ 
tage  of  fire  modeling  is  the  significantly  reduced  need  for  full-scale  fire 
tests  which  are  extremely  expensive  and  time  consuming. 

Fire  models  can  be  categorized  as  either  zone  models  or  field  models.  Zone 
models  deal  with  dividing  the  fire-affected  environment  or  space  into  distinct 
regions  or  zones  which  can  be  separately  analyzed  either  empirically  or  theore¬ 
tically  in  terms  of  input  and  output  information  based  on  mass  and  energy  balan¬ 
ces.  Examples  of  such  zones  include  the  fire  envelope,  fire  plume,  hot  gas 
layer,  ceiling  jets,  wall  jet,  flow  through  openings,  etc.  These  input  and 
output  quantities  from  individual  zones  are  then  assembled  to  describe  the 
overall  gas  dynamic  phenomenon  of  the  entire  fire  environment.  In  most  cases, 
there  results  a  set  of  ordinary  differential  equations  in  time  which  are  then 
solved  numerically  by  standard  integration  algorithms.  While  the  zone  models 
are  generally  computationally  efficient,  shortcomings  do  exist  in  that  the 
models  for  some  of  the  zones  are  not  adequately  known  and  quantified,  and  some¬ 
times  there  is  a  question  as  to  exact  range  of  validity  of  such  zone  models. 
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Efforts  are  now  being  continued  to  refine  the  existing  zone  models  and  also  to 
compare  with  full-scale  test  data  for  validation  purposes  (for  example  see 
Rockett  et  al.,  1987). 

Field  models  which  are  based  on  the  numerical  solutions  to  the  governing 
differential  field  equations  of  the  conservation  of  momentum,  mass,  energy  and 
spices,  on  the  other  hand,  have  not  been  widely  used  in  the  past  for  fire 
modeling  purposes  even  though  these  models  are  inherently  more  rational  and 
capable  of  revealing  both  the  important  large-and  small-scale  phenomena  in  a 
given  fire-spread  problem.  The  primary  reasons  are  the  excessive  use  of  com¬ 
puting  resources  and  degrees  of  uncertainty  in  modeling  such  phenomena  as  tur¬ 
bulence,  combustion  and  thermal  radiation.  However,  particularly  because  of  the 
increasing  availability  of  high-performance  computers  including  supercomputers, 
field  models  for  fire  modeling  have  found  increasing  use  in  more  recent  times. 
For  instance,  such  models  have  been  developed  at  the  National  Bureau  of 
Standards  (Baum  and  Rehm,  1984),  Imperial  College  of  London  (Bagnaro,  Laouisset 
and  Lockwood,  1983),  Borehamwood  Fire  Research  Station  (Markatos  et  al.,  1982), 
and  the  University  of  Notre  Oame  (Yang  et  al.,  1984,  Yang  and  Lloyd,  1985),  Kou 
et  al.,  1986,  and  Nies,  1986).  In  view  of  the  different  physical  and  chemical 
models  utilized  in  each  of  these  field  models  and  the  fact  that  differences  also 
exist  in  the  details  of  the  numerical  algorithms  involved,  there  is  no  way  to 
ascertain  the  relative  accuracies  of  these  models.  Comparisons  of  the  numerical 
results  with  full-burn  tests  are  scarce,  but  do  exist  and  are  fairly  reasonable 
(Yang  et  al.,  1984;  Yang  and  Lloyd,  1985;  Bagnaro  et  al.,  1983,  Markatos  and 
Pericleous,  1983).  An  ultimate  assessment  of  the  various  field  models  cannot  be 
made  until  more  validation  studies  based  on  full-scale  experimental  data  become 
available.  Additional  discussions  on  fire  modeling  have  been  recently  given  by 
Stroup  (1987). 
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Despite  the  above-mentioned  uncertainties  of  the  field  models,  many  of  the 
important  physical  effects  can  now  be  accommodated.  Those  include  turbulence 
(Vang  and  Liu,  1978;  Bagnaro  et  al.,  1983;  Marketos  and  Pericleous,  1983),  ther¬ 
mal  radiation  (Lloyd  et  al.,  1979;  Bagnaro  et  al.,  1983;  Markatos  and 
Pericleous,  1983);  Nies,  1986),  combustion  (Bagnaro  et  al.,  1983;  Markatos  and 
Pericleous,  1983),  three-dimensional  vented  spaces  (Bagnaro  et  al.,  1983; 

Markatos  and  Pericleous,  1983;  Yang  and  Lloyd,  1985;  Nies,  1986),  closed-space  in¬ 
duced  pressure  rises  (Nies,  1986),  and  wall  losses  (Nies,  1986).  However,  all  the 
above-mentioned  field  models  deal  with  box-like  rectangular  spaces  or  enclo¬ 
sures,  and  yet,  many  full-burn  facilities  have  other  geometries  for  which  these 
field  models  are  very  difficult  to  apply.  A  good  example  is  the  facility  known 
as  FIRE  I  located  at  the  Naval  Research  Laboratory  (NRL) ,  which  is  in  the  form 
of  a  pressure  vessel  including  a  central  cylindrical  body  plus  two  spherical  end 
caps  (Alexander  et  al.,  1982).  The  purpose  of  this  report  is  to  present  another 
field  model  which  is  based  on  a  generalized  orthogonal  coordinate  system  which 
still  retains  the  same  capabilities  in  terms  of  physical  models  for  turbulence, 
full  compressibility,  strong  buoyancy,  surface-to-surface  and  f lame-to-surface 
radiation,  wall  losses,  and  effects  of  heat  addition  in  a  closed  space.  The  rec¬ 
tangular  coordinate-system  case  as  given  by  Nies  (1986)  becomes  a  special  case 
of  this  new  field  model,  and  it  can  be  mentioned  that  the  model  of  Nies  (1986) 
is  basically  an  extension  of  the  original  two-dimensional  UNDSAFE-II  code  devel¬ 
oped  at  the  University  of  Notre  Dame  (Yang  and  Liu,  1978)  for  vented  rec¬ 
tangular  enclosures.  Furthermore,  this  new  field  model  will  be  demonstrated  by 
presenting  a  set  of  numerical  results  which  simulate  some  experimental  data 
obtained  in  the  NRL  FIRE  I  facility.  Finally,  several  additional  capabilities 
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of  this  new  field  model  designed  to  better  simulate  reel  fire  scenarios  will 
also  be  delineated. 


MATHEMATICAL  FORMULATION  OF  THE  FIELD  MODEL 

As  pointed  out  previously,  all  field  models  for  the  fire-spread  phenomena 
are  based  on  first-principle  conservation  differential  equations  for  mass, 
momentum,  energy  and  species.  The  approach  utilized  here  is  similar  to  that 
given  by  Yang  et  al.  ( 1987a, b),  and  basically,  the  governing  equations  in  the 
Cartesian  coordinates  are  first  transformed  into  those  of  a  generalized  orthogo¬ 
nal  coordinate  system,  which  are  then  f inite-differenced  by  the  control -volume 
staggered-cel 1  primitive-variable  approach.  The  flux  terms  are  handled  in  such 
a  way  that  the  final  finite-difference  algebraic  equations  resemble  closely 
those  of  the  Cartesian  coordinates.  As  a  result,  the  computer  algorithm  for  the 


generalized  orthogonal  coordinates  is  no  more  complex  than  that  for  the 
Cartesian  coordinates.  In  fact,  it  has  recently  been  shown  that  the  present 
approach  represents  a  viable  alternative  to  the  popular  approach  of  utilizing 
body-fitted  coordinates  (Yang  et  al.,  1988).  Under  the  conditions  of  turbulent 
compressible  flow  with  buoyancy,  the  governing  equations  for  the  mean  quantities 
in  the  Cartesian  coordinates  may  be  written  as  follows  in  the  standard  tensor 
forms: 


°t  ♦  (PUi),! 


(pui)t  +  (PUiUjJ^j  -  -p#i  -  pGi  + 


(PCpaT)t  +  (PUiCpnT)^  —  qciJ  +  Q 


■ 
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Whert  smoke  propagation  in  the  space  needs  to  be  considered,  an  additional 


equation  similar  to  (3)  for  the  smoke  concentration  can  be  added  (Yang  and  Liu, 


1978;  Raycraft,  1987).  It  is,  however,  not  shown  here,  since  the  smoke  field  is 


not  considered  in  the  demonstration  of  this  field  model,  as  will  be  shown 


later,  and  in  the  comparison  of  the  numerical  results  and  the  test  data.  In 


Eqs.  (1),  (2)  and  (3),  the  physical  coordinates  are  with  i  =  1,  2  and  3, 
p  is  the  mean  fluid  density,  u..  are  the  mean  velocity  components  in  x^, 


subscript  t  signifies  the  time  derivatives,  p  is  the  mean  static  pressure. 


represents  the  gravity  vector  components,  T  is  the  mean  fluid  temperature, 


and  Q  represents  the  source  term.  It  is  noted  that  inside  the  fire  envelope. 


Q  includes  both  the  thermal  radiation  and  combustion  heat  contributions,  and 


outside  the  envelope,  Q  simply  vanishes,  provided  that  the  gas  radiation 


effects  are  neglected.  Furthermore,  the  mean  isobaric  heat  capacity  c  ,  the 

r'" 


shear  stress  tensor  o^,  and  the  conduction-flux  qci  are  defined  as 


'pm  ~  T 


I 

"-TT  I  cp(WT 


Oij  -  ueff  (Ui,j  +  Uj,i  -  2/3 


qci  =  -keff  T,i 


respectively,  where  TR  is  a  convenient  reference  temperature  generally  taken 


to  be  the  fluid  temperature  prior  to  the  commencement  of  the  fire,  5^  is  the 
usual  Kronecker  delta,  and  neff  and  ke^  are  the  effective  dynamic  viscosity 
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and  thermal  conductivity,  respectively,  both  accounting  for  the  combined  lami¬ 
nar  and  turbulent  contributions.  As  will  be  shown  later,  both  pg^  and 
need  to  be  modeled  separately. 

In  addition,  the  fluid  medium  is  taken  to  be  a  perfect  gas,  and  the  perfect 
gas  relation 

P  =  pRT  (7) 

where  R  is  the  fluid  gas  constant,  is  not  only  used  to  accommodate  the 
compressibility  effects,  but  also  needed  to  account  for  pressure  build-up  due  to 
heating  from  the  fire  in  a  constant-volume  space.  As  an  approximation,  all 
laminar  transport  properties, such  as  cpm,  p  and  k  of  the  fluid, are  considered 
to  be  temperature  independent  and  calculated  at  TR.  This  is  justifiable  because 
c  "is  normally  only  a  weak  function  of  temperature,  while  p  and  k  only 

pin 

contribute  to  a  minor  degree  to  the  respective  effective  values  due  to  the  much 
larger  turbulent  effects  expected  in  the  fire  phenomena.  The  initial  conditions 
are  simply  that  there  is  no  motion  of  the  fluid,  i.e.,  u.  =  0,  and  T  =  T  . 

1  K 

The  boundary  conditions  are  the  usual  no-slip  conditions  at  the  enclosure 
interior  wall  surfaces  and  that  for  the  temperature  field  is  determined  by  a 
heat  balance  there.  When  the  interior  wall  is  heated  up,  heat  loss  occurs 
through  the  wall  thickness  and  eventually  to  the  environment  from  the  exterior 
wall.  In  the  present  study,  such  wall  losses  are  accommmodated  by  transient  con¬ 
duction  through  the  wall  and  then  by  combined  convection  and  radiation  at  the 
outer  surface  by  means  of  a  constant  coefficient  of  heat  transfer.  Since  the 
wall  thickness  is  always  small  compared  to  the  characteristic  length  along  the 
wall,  a  one-dimensional  conduction  analysis  should  be  sufficient.  Thus,  for 
conduction  through  the  wall  thickness,  we  may  write 
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ps  cs  Tst  =  (MJ 
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where  the  subscript  refers  to  the  wall  material,  and  n  is  the  outward  normal 
coordinate,  together  with  the  boundary  conditions: 


n  =  0 


n  =  L 


n  1/  ^  I  i  S 

qr  "  keff  an  =  'ks  aTT 


‘ks  3?T  *  h<T  -  Te) 


where  n  =  0  is  at  the  inner  surface,  Lw  is  the  wall  thickness,  h  is  the 
coefficient  of  heat  transfer  at  the  outer  surface,  Te  is  the  environment  tem¬ 
peratures  and  qp  is  the  thermal  radiation  flux  arriving  at  the  surface  from 
the  enclosure  interior.  As  pointed  out  previously,  h  is  taken  to  be  a  con¬ 
stant,  depending  on  the  wind  conditions  outside.  Before  Eqs.  (1),  (2),  (3) 
and  (8)  can  be  solved  simultaneously,  additional  turbulence  and  thermal  radia¬ 
tion  models  must  be  introduced,  and  these  models  will  be  given  in  a  later  sec¬ 
tion. 

Now  the  above  governing  equations  in  x..  are  ready  to  be  transferred 
into  those  of  a  generalized  orthogonal  coordinate  system  given  by  the  coor¬ 
dinates  8 1 .  The  two  coordinate  systems  are  related  by  means  of  the  scale 
factors  h^  in  the  direction  01  given  by 


hi  =  (g-gi) 


ax.  ax.  1/2 


(aex  ael) 
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where  is  the  base  vector.  One  of  the  special  properties  of  the  orthogonal 
coordinates  is  that  these  base  vectors  are  orthogonal  so  that  the  covariant  and 
contravariant  metric  tensors  satisfy  the  following  conditions 

«ij  -  5i  •  5j  -  Sijhihj  (12) 


’i3-o£ 


resulting  in  a  diagonalized  metric  tensor,  where  g  is  the  determinant  of  g^., 


9  •  l9ij I  »  h|  hj  h§ 


With  the  above  orthogonal  properties,  transformations  of  Eqs.  (1),  (2)  and 
(3)  in  the  x..  system  lead  to  the  following  respective  equations  in  the  91 
system: 


pt  +  T7T  rr(9l/2  p  ul/V  *  0 

g  do 


(0“l)t  *  ZI75  uV/V  "  '  p.i/hi  +  +  -rn  ^rigVJcK 
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where 
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The  above  conservation  equations  are  clearly  valid  for  any  orthogonal  coor¬ 
dinate  systems.  In  particular,  h1  =  h2  =  h3  =  1  for  the  Cartesian  coor¬ 
dinates;  h1  =  r  =  0*,  h2  =  h3  =  1  for  the  cylindrical  coordinate  system  where 
r  is  the  radial  variable;  and  h1  =  rsinij)  =  0*sin93,  h2  =  1,  h3  =  r  =  0*  for 
the  spherical  coordinates,  where  $  is  the  polar  angle.  The  boundary  and  ini¬ 
tial  conditions  need  no  further  transformation.  In  addition,  in  view  of  the 
relatively  thin  wall  thickness,  one-dimensional  conduction  is  still  valid  and 
therefore,  Eqs.  (8),  (9)  and  (10)  remain  the  same. 

As  pointed  out  previously,  before  the  governing  equations  (15),  (16),  (17) 
and  (8)  can  be  solved,  physical  models  for  the  radiation  exchange  in  the  enclo¬ 
sure  and  for  the  turbulence  field  are  still  needed.  Since  the  radiation 
exchange  depends  on  the  specific  geometry  of  the  enclosed  space,  the  radiation 
model  will  be  described  later  after  the  geometry  of  the  NRL  FIRE  I  test  facility 


is  introduced.  However,  a  turbulence  model  utilized  in  the  present  study  can  be 
given  here,  again  in  terms  of  an  arbitrary  orthogonal  coordinate  system.  This 
model  is  based  on  a  mixing-length  formulation  originally  developed  for  two- 
dimensional  recirculating  flows  (Nee  and  Liu,  1978),  and  has  been  extended  to 
three-dimensional  flows  with  good  success  (Yang  and  Lloyd,  1985).  In  terms  of 
the  orthogonal  coordinate  system,  the  equation  for  the  effective  viscosity 
according  to  this  model  can  be  written  as 
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where  pR  is  the  laminar  viscosity  at  TR,  Pr^  is  the  turbulent  Prandtl 
number  taken  to  be  unity,  and  R,.  is  the  gradient  Richardson  number  given  by 
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where  UR  is  a  reference  veloci  ty,  and  n  is  a  unit  vector  in  the  negative 
gravity  direction.  In  addition,  the  mixing  length  I  normalized  with  the 
enclosure  height  H  can  be  written  as: 
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where  K  is  a  constant,  normally  taken  to  be  0.2.  Finally,  the  effective 
conductivity  ke^ 


is  given  by 
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where  kR  is  the  molecular  thermal  conductivity  at  TR,  and  Pr  is  the 
corresponding  laminar  Prandtl  number. 


NUMERICAL  FINITE-DIFFERENCE  CALCULATIONS 

In  the  present  study,  the  governing  equations  (15),  (16)  and  (17)  in  the 
generalized  orthogonal  coordinates  are  solved  by  means  of  the  primitive-variable 
staggered  control -volume  finite-difference  approach.  As  recently  pointed  out  by 
Vang  et  al.  (1987b,  1988),  one  major  advantage  of  utilizing  the  fully  transformed 
equations  is  that  by  judiciously  choosing  new  definitions  of  flux  and  source 
quantities,  the  final  finite-difference  equations  appear  very  much  like  those 
for  the  well-known  Cartesian-coordinate  cases,  and  hence  can  be  handled  in  a 
routine  manner  for  their  solutions.  That  this  is  the  case  can  be  illustrated  as 
follows: 

For  a  control  volume  given  by  AV  =  /g  A01AS*A03,  where  A  indicates  a 
step  size,  the  mom  entum  equation  (16)  for  u1,  after  integrated  over  the 
staggered  control  volume  surrounding  a  grid  point  P,  becomes 
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and  M1J  is  the  total  momentum  flux  along  the  0J  direction  for  the  velocity 
component  u1  due  to  convection  and  diffusion.  The  subscripts  e,  w,  n,  s,  f 
and  b  denote  the  east,  west,  north,  south,  front  and  back  boundaries  one-half 
cell  size  away  from  the  point  P,  respectively.  More  specifically, 
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Here  the  idea  of  Raithby  et  al.  (1986)  of  using  a  stress-flux  formulation  is 
utilized,  and  the  momentum  flux  is  then  split  into 
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As  a  result,  Eq.  (23)  now  becomes 
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Here  the  stresses  are  evaluated  from  the  information  of  the  prior  itera¬ 


tion,  and  the  source  term  is  known  at  the  present  iteration.  Once  these  are 
properly  evaluated,  Eq.  (30)  becomes  identical  in  form  to  that  for  the 


Cartesian  coordinates.  Further  approximations  to  M1J  can  be  made.  Basically, 
the  convective  terms  are  approximated  by  the  QUICK  scheme  of  Leonard  (1983), 
extended  to  three-dimensional  cases,  and  the  diffusion  terms  are  by  central  dif¬ 
ferences.  The  details  have  been  given  by  Raycraft  (1987),  and  hence  will  not  be 
repeated  here. 

The  treatment  of  the  energy  equation  is  very  similar.  For  example,  a  total 
heat  flux  J1,  which  includes  both  convection  and  conduction,  is  introduced  as 


,1 
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The  corresponding  finite-difference  equation  for  the  energy  equation  (17)  then 
becomes 
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(p  Cpm  T)^  AV  +  Je^Ae  -  Jw'aw  +  Jn2An 

-  Js2As  ♦  Jf3Af  -  Jb3Afa  =  ST  (33) 

where  the  source  S^.  also  includes  Q  and  the  solution  procedure  is  the  same 
as  that  for  the  momentum  equations. 

In  addition,  it  is  noted  that  for  a  closed  enclosure,  any  energy  addition 
to  the  gas  due  to  the  fire  tends  to  raise  the  pressure  in  the  enclosure  in  view 
of  the  constant  volume.  Consequently,  a  time-dependent  global  pressure 
correction  is  needed  to  determine  this  rise,  in  addition  to  the  local  pressure 
correction  in  the  primitive-variable  control -volume  procedure,  which  only 
depends  on  the  local  velocity  variations.  For  this  purpose,  the  technique  based 
on  the  perfect  gas  relation  as  developed  by  Nicolette  et  al.  (1985)  is  utilized 
in  the  present  study,  and  this  technique  has  also  been  described  by  Raycraft 
(1987). 

NRL  FIRE  I  TEST  FACILITY  AND  RADIATION  MODEL 

As  pointed  out  previously,  the  validation  of  the  field  model  developed  in 
the  present  study  is  based  on  test  data  obtained  in  the  Naval  Research 
Laboratory  (NRL)  FIRE  I  full-scale  test  facility  (Alexander  et  al.,  1982).  As 
shown  schematically  in  Fig.  1,  the  facility  has  a  cylindrical  midsection  with  a 
radius  of  9.6  feet  (2.93m)  and  a  length  of  27.4  feet  (8.352m)  and  hemispherical 
end  caps.  Thus,  it  is  in  the  form  of  a  pressure  vessel  having  a  total  volume  of 
11,639  cubic  feet  (324m3).  The  instrumentation  consists  of  pressure  trans¬ 
ducers,  thermocouples,  and  radiometers.  There  is  also  capability  to  measure 
smoke  obscuration  levels,  gas  composition  and  humidity.  The  use  of  circulation 
fans  can  also  be  included  in  a  test  to  determine  the  effect  of  ventilation. 
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Pressure  transducers  are  located  at  the  north  and  south  ends  of  the  test 
chamber,  while  the  temperatures  are  monitored  with  thermocouple  arrays  located 
inside  the  spherical  end  caps  as  shown  in  Fig.  1.  These  are  chromel-alumel 
thermocouples  with  diameters  of  0.12  mm  and  have  ceramic  insulation  enclosed  in 
304  stainless  steel  jackets  1.0  mm  in  diameter.  As  will  be  shown  later,  these 
temperature  readings  provide  the  data  for  comparison  with  the  numerical  calcula¬ 
tions  based  on  the  field  model. 

In  the  field  of  fire  modeling,  the  importance  of  thermal  radiation  contri¬ 
bution  to  the  energy  transfer  and  turbulent  buoyant  flow  processes  involved  is 
generally  well  recognized.  However,  accommodation  of  the  radiation  effects  in 
fire  models  is  always  difficult  and  requires  solutions  to  the  radiative 
transfer  equation  for  complex  geometries  and  nongray-gas,  nonhomogeneous ,  and 
nonisothermal  behaviors  (Yang,  1986).  The  presence  of  soot  and  smoke  further 
complicates  the  situation,  even  though  some  simplification  here  is  possible. 
Therefore,  any  reasonable  simplification  in  including  the  radiation  effects  is 
desirable.  In  the  present  study,  only  surface-surface  and  flame-surface 
radiation  exchange  is  included  in  the  field  model.  In  other  words,  the  gas  out¬ 
side  the  flame  envelope  is  taken  to  be  transparent.  In  addition,  the  flame  is 
considered  to  be  represented  by  an  equivalent  gray  surface  which  exchanges 
radiation  with  all  other  surfaces  in  the  enclosure.  This  simplification  is 
deemed  reasonable,  since  it  is  known  that  the  effects  of  participating  gases  in 
an  enclosure  is  largely  to  equalize  the  gas  temperatures,  and  the  radiation 
exchange  is  in  general  dominated  by  surface-surface  interactions  (Chang  et  al., 
1983).  Under  this  simplification,  the  radiation  heat  transfer  can  be  readily 
calculated  by  any  one  of  the  well-known  enclosure  procedures  (Sparrow  and 


Cess,  1978;  Siegel  and  Howell,  1981).  In  the  present  study,  the  radiation 
calculations  are  based  on  the  net  radiosity  method  which  is  particularly  suited 


for  gray  surfaces.  The  radiosity  B^  at  a  surface  i  can  be  written  as 

8i  *  ciaTj  +  0  -  e.)  I  B.F. .  <3<) 

’  '  ’  j-1  J  ’J 

where  e  is  the  emissivity,  and  F.^  is  the  direct  view  factor  from  surface  i 
to  surface  j.  The  index  N  denotes  the  total  number  of  surfaces  in  the  enclo¬ 
sure,  including  that  of  the  flame.  The  above  equation  results  in  N  linear, 
inhomogeneous,  algebraic  equations  for  the  N  unknown  radiosities.  By  solving 
the  simultaneous  equations,  usually  by  matrix  inversions,  B .  can  be  found, 
and  the  radiation  heat  flux  at  surface  i  can  then  be  readily  calculated  from 


(35) 


which,  for  all  enclosure  interior  surface  elements,  becomes  part  of  the  thermal 
boundary  conditions  in  accordance  with  Eq.  (9).  This  method  of  calculating  q^.. 
can  be  further  simplified  by  obtaining  a  direct  relation  between  qpi  and  the 
temperatures  T..  so  that  matrix  manipulations  can  be  kept  at  a  minimum.  Eq. 
(34)  can  be  rewritten  as 
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where 


which,  in  a  matrix  form,  becomes 


[X]  <B>  -  a<T4> 


(38) 


By  inverting  the  coefficient  matrix  [X],  we  obtain  for  the  radiosity 
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Then  in  accordance  with  Eq.  (35),  we  finally  have 
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Here  X^-1  only  depends  on  the  emissivities  which  are  taken  to  be  constants, 
and  consequently,  together  with  G^j,  can  be  calculated  once  for  all  in  a  given 
problem. 

It  is  thus  clearly  seen  that  the  radiation  model  utilized  in  the  present 
study  essentially  reduces  to  a  problem  of  determining  all  the  appropriate  view 
factors  F^j.  To  be  consistent  to  the  finite-difference  calculations  for  the 
buoyant  flow  in  the  enclosure,  in  which  each  interior  surface  cell  has  a  tem¬ 
perature  which  depends  on  time,  these  surface  cells  can  also  be  conveniently 
used  as  surfaces  i  for  radiation.  Consequently,  the  determination  of  the  view 
factors  does  not  require  any  surface  integration,  and  they  can  be  calculated 
directly  from  the  definition 
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where  0  is  the  angle  between  the  normal  of  the  surface  and  R,  the  line  con¬ 
necting  and  Aj.  Since  depends  strongly  on  the  geometry  of  the 

enclosure,  no  generic  forms  can  be  written  in  terms  of  the  generalized  ortho- 
gonala  coordinates,  and  this  is  the  very  reason  that  the  radiation  model  cannot 
be  introduced  until  the  geometry  is  specified. 

For  the  geometry  of  FIRE  I,  the  view  factors  F^  depend  on  such  surface 
cell  pairs  as  internal  cylindrical  to  cylindrical,  cylindrical  to  hemispherical, 
and  hemispherical  to  hemispherical  surfaces.  For  the  last  case,  surfaces  i 
and  j  may  also  be  located  on  different  end  caps.  Surprisingly,  not  all  such 
view  factors  are  known  in  the  literature  (Siegel  and  Howell,  1981).  As  a 
result,  these  view  factors  have  been  derived  directly  from  the  FIRE  I  geometry. 
They  have  been  given  by  Raycraft  (1987)  and  there  is  no  need  to  repeat  them 
here.  Furthermore,  the  view  factors  associated  with  the  flame  must  be  treated 
somewhat  differently,  and  this  will  be  given  in  the  following  section. 


RESULTS  OF  NUMERICAL  CALCULATIONS  AND  DISCUSSIONS 

Detailed  numerical  calculations  based  on  the  field  model  have  been  carried 
out  for  the  NRL  FIRE  I  facility  geometry  so  that  the  results  can  be  compared 
directly  with  the  test  data.  One  set  of  test  data  for  burning  methanol  inside 
the  FIRE  I  tank  is  chosen  for  this  purpose.  For  this  burn  test,  the  tank  is 
essentially  empty  and  the  fuel  pan  is  located  at  23.1  feet  (7.041m)  from  either 
end  cap  and  elevated  3.21  feet  (0.978m)  from  the  bottom.  The  initial  tem¬ 
perature  and  pressure  are  at  35.6*C  and  1  atm.,  respectively.  Shown  in  Fig.  2 
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is  the  grid  system  utilized  in  the  computations.  As  seen,  there  are  20 
segments  in  the  circumferential  direction,  12  radial  segments  representing  the 
intorior  of  the  tank,  and  18  segments  along  the  cylindrical  midsection,  for  a 
total  of  6,720  interior  cells.  Also,  in  order  to  avoid  calculation  difficulties 
surrounding  the  line  of  zero  radius,  a  row  of  cells  is  added  to  preserve 
the  direct  of  the  velocity  vectors  (Yang  et  al . ,  1987b).  The  thickness  of 
the  wall  is  0.375  in  (0.95  cm),  and  occupies  one  single  cell  segment. 
Consequently,  the  entire  interior  surface  of  the  tank  has  560  surface  cells 
which  all  exchange  radiation  with  the  flame  and  also  with  each  other,  and  there 
are  also  560  cells  for  conduction  through  tie  wall  thickness.  Furthermore, 
simi .ar  to  chat  in  the  UN0SAFE-I1  field  mocel  (Yang  and  Liu,  1978),  the  flame  is 
modeled  as  *  volumetric  heat  source,  and  in  the  present  calculations,  the  entire 
heat  release  rate  in  the  flame  is  taken  to  be  distributed  uniformly  over  19 
cells  above  the  fuel  pan.  Thus,  radiation  heat  transfer  calcula-tions  involve 
altogether  579  surfaces.  The  numerical  calculations  also  require  several  physi¬ 
cal  data  relative  to  the  wall  properties  and  external  boundary  conditions.  In 
addition  to  the  wall  thickness  given  previously,  the  specific  heat,  thermal  con¬ 
ductivity,  and  density  of  the  metallic  wall  (ASTM-285  Grade  C  steel)  are  taken 
to  be  0.1  B/( lbF)  (0.42  kJ/kgK),  25  B/hr.ft.F(43.27  W/mK),  and  487  lb/ft» 

(7,801  kg/a*),  respectively.  There  is  also  convection  and  radiation  heat 
transfer  occurring  at  the  tank  exterior  surface,  and  the  corresponding  coef¬ 
ficient  of  heat  transfer  is  taken  to  be  a  constant  of  15  B/hr.ft*F  (85.169 
W/m*K).  This  coefficient  can  be  adjusted  in  the  computations,  if  needed.  The 
time  increment  is  essentially  determined  by  numerical  stability  requirement. 
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The  driving  force  in  the  physical  problem  is  the  heat  release  rate  from  the 
fire.  Since  in  the  present  model  the  flame  is  taken  to  be  a  volumetric  heat 
source,  the  unsteady  heat  release  rate  must  be  treated  as  an  input.  Normally, 
such  data  can  be  obtained  from  the  fuel  depletion  rate  and  the  heating  value  of 
the  fuel.  Unfortunately,  because  of  instrumentation  error  during  the  methanol 
fire  test,  the  fuel  depletion  rate  data  were  also  in  error,  and  consequently  the 
heat  release  rate  must  be  determined  in  an  indirect  way.  Since  the  static 
pressure  in  the  tank  during  the  fire  test  was  recorded,  and  in  a  constant  volume 
heat  addition  process,  the  heat  added  is  directly  proportional  to  the  rate  of 
rising  static  pressure,  an  estimate  of  the  heat  release  rate  can  be  obtained 
from  the  recorded  pressure  data.  In  an  attempt  to  generate  the  unsteady  heat 
release  rate  in  this  way  automatically  on  the  computer,  a  pressure-tracking 
routine  based  on  feedback  control  and  the  field  model  has  been  developed 
(Nies,  1986;  Raycraft,  1987).  The  result  is  shown  in  Fig.  3.  Even  though  the 
routine  is  performing  well  insofar  as  the  pressure  level  is  concerned,  the 
resulting  heat  release  rate,  which  is  sensitive  to  the  time  derivative  of  the 
pressure,  oscillates  with  large  variations  leading  also  to  large  variations  in 
the  temperature  field  in  the  tank.  This  heat  release  rate  from  the  pressure¬ 
tracking  is  given  in  Fig.  4.  Since  such  large  variations  of  the  heat  release 
rate  are  not  expected  in  the  fire  test,  it  is  reasonable  to  use  a  curve-fitted 
variation  as  the  input  to  the  field  model  calculations,  as  also  shown  in  Fig.  4. 
In  order  to  demonstrate  that  such  a  curve-fitted  heat  release  rate  variation  is 
a  reasonable  approximation  to  the  real  data,  the  calculated  pressure  rise  is 
also  given  in  Fig.  3,  and  it  is  seen  that  it  does  follow  the  experimental  curve 
rather  well.  (Detailed  intermediate  results  based  on  pressure  tracking  can  be 
found  in  the  thesis  of  Raycraft  (1987). 
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Figure  3.  Pressure  History  After  Commencement  of  Fire 


HEAT  RELEASE  RATE  (  KW ) 


A  more  important  result  of  the  numerical  computations  is  to  compare  the 
temperatures  at  the  various  thermocouple  locations  with  those  from  the  test. 


The  thermocouples  are  cl  1  mounted  on  vertical  racks  in  the  two  end-cap  regions 
close  to  the  edges  of  the  cylindrical  midsection  (Fig.  1).  The  three  ther¬ 
mocouples,  1,  2  and  4,  are  located  in  the  south  hemispherical  end  cap,  and  the 
readings  there  are  compared  with  the  numerical  results.  Thermocouple  1  is 
located  79  inches  (2.01  m)  above  the  midplane  of  the  tank.  Thermocouple  2  is  one 
foot  (0.305  m)  below  thermocouple  1,  and  thermocouple  4  is  at  two  feet  (0.61  m) 
below  thermocouple  2.  The  thermocouple  3  was  defective  during  the  test,  and 
hence  comparisons  there  are  not  made.  Figs.  5,  6  and  7  show  the  temperature 
variations  at  thermocouples  1,  2  and  4  locations,  respectively,  as  compared  to 
those  from  the  test,  which  represent  smoothed  mean  temperatures.  The  agreement 
is  reasonable  at  the  two  lower  thermocouple  locations  with  maximum  deviations  of 
less  than  10#C.  However,  the  calculations  definitely  overpredict  the  tem¬ 
peratures  at  the  thermocouple  1  location  by  a  maximum  of  about  30*C.  Several 
possible  reasons  may  be  responsible.  One  is  that  since  the  thermocouple  1  loca¬ 
tion  is  closest  to  the  tank  wall,  its  temperature  is  most  affected  by  the  boun¬ 
dary  condition  at  the  tank  exterior  surface  in  terms  of  outside  ambient 
temperature  and  the  coefficient  of  heat  transfer,  both  of  which  have  been  estimated  for 
calculation  purposes.  Another  possible  reason  is  that  in  the  model,  the  heat 
release  rate  is  taken  to  be  uniformly  distributed.  In  reality,  considerable 
vertical  nonuniformity  may  result  from  the  combustion  process  in  the  flame.  If 
more  heat  is  released  in  the  lower  part  of  the  flame,  more  heat  would  be 
radiated  away  there  so  that  the  thermal  plume  would  attain  a  lower  temperature, 
resulting  also  in  lower  temperatures  at  the  thermocouple  1  location.  A  third 


Thermocouple  #1 


possible  reason  is  that  gas  and  soot  radiation  is  neglected  in  the  field  model. 
Since  in  general  the  effect  of  gas  and  soot  radiation  is  to  make  the  tempera¬ 
ture  field  more  uniform  in  the  enclosure  (Yang,  1986),  the  thermocouple  1  loca¬ 
tion  would  also  tend  to  have  lowet  temperatures.  Also,  the  uncertainty  in  the 
experimental  heat  release  rate  itself  could  also  be  responsible  for  the  degree 
of  discrepancy  shown  in  Fig.  5.  Incidentally,  it  is  noted  in  Figs.  5,  6  and  7 
that  all  calculated  temperatures  exhibit  slow  oscillations  which  cannot  be 
attributed  to  turbulence.  Such  oscillations  could  very  well  be  real  because 
of  inherent  thermal  instability  associated  with  thermal  plumes  (Yang  et  al., 
1983). 

One  of  the  major  advantages  of  the  field  model  is  the  fact  that  the  numeri¬ 
cal  solution  gives  detailed  spatial  and  temporal  variations  of  the  velocity 
and  temperature  fields.  Such  information  is  critical  in  the  development  of 
fire  mitigation  measures,  especially  in  the  placement  of  barriers  and  par¬ 
titions,  and  in  the  determination  of  fire  loads  on  the  walls.  The  field 
results  from  the  calculations  in  terms  of  the  velocities  and  temperaturess  can 
now  be  shown.  Because  of  the  three-dimensional  nature  of  the  fields,  it  is 
necessary  to  look  at  these  quantities  at  specific  sections  of  the  tank.  At  a 
time  instant  30  seconds  into  the  fire,  the  isotherms,  or  the  temperature  field, 
and  the  velocity  vector  field  are  shown  in  Figs.  8.  Fig.  8(a)  gives  a  front 
view  of  these  fields  at  the  mid  section  of  the  tank.  Fig.  8(b)  refers  to  the 
view  also  at  the  mid  section  from  the  south  end  of  the  tank,  while  Fig.  8(c) 
is  another  view  from  the  south  end,  but  at  the  intersection  between  the  end 
cap  and  the  cylindrical  mid  section,  or  at  the  base  of  the  end  cap.  Some 
sections  are  used  in  Figs.  9,  10,  11  and  12  at  60  seconds,  90  seconds,  120 
seconds  and  150  seconds  from  the  commencement  of  the  fire,  respectively. 


Section  View  at  Base 


Since  the  fire  is  located  at  the  center  of  the  tank,  both  sections  in 
Figs.  8{a)  and  8(b)  include  the  flame,  and  therefore  it  is  seen  that  the  fire 
plume  is  already  well  formed  at  30  seconds.  There  is  also  a  hot  ceiling  layer 
along  the  ceiling  of  the  tank,  and  the  lower  two  thirds  of  the  tank  are  still 
at  almost  the  initial  temperature.  The  general  flow  circulation  pattern  is 
much  more  complex  than  might  be  expected.  The  entrainment  into  the  fire  plume 
is  limited  to  the  immediate  neighborhood  of  the  fire,  as  can  be  seen  in  Figs. 
8(a)  and  8(b).  The  flow  in  the  hot  ceiling  layer  does  not  seem  to  have  strong 
enough  momentum  to  carry  the  flow  toward  the  lower  half  of  the  tank,  but  re¬ 
turns  into  the  tank  interior,  resulting  in  a  downward-biased  flow  even  at  the 
section  as  shown  in  Fig.  8(c).  The  flow  completes  its  journey  back  toward  the 
fire  region  in  a  somewhat  spatially  oscillatory  path.  At  60  seconds  as  shown 
in  Fig.  9,  the  depth  of  the  hot  ceiling  layer  is  increased,  even  though  the 
lower  half  of  the  tank  is  still  not  affected  much  by  the  fire.  The  flow  in 
the  ceiling  layer  toward  and  down  the  spherical  end  cap  is  definitely  spiral 
in  nature.  In  addition,  the  spatially  oscillatory  flow  back  to  the  fire  zone 
become  sufficiently  large  that  it  induces  a  backward  flow  toward  the  end  cap 
again  at  the  bottom  of  the  fire.  It  may  also  be  of  interest  to  note  that  the 
isotherms  are  now  more  densely  packed  at  the  ceiling  region  of  the  tank,  indi¬ 
cating  that  heat  losses  through  the  wall  now  become  more  important.  Similar 
trends  of  the  temperature  and  velocity  fields  can  be  seen  at  90,  120  and  150 
second  instants  in  Figs.  10,  11  and  12,  respectively.  It  is  interesting  to 
note  that  the  hot  ceiling  layer  reaches  down  almost  to  the  half  tank  height 
level  at  120  and  150  seconds,  and  inside  the  hot  layer,  the  temperature  field 
is  essentially  stably  stratified  except  a  small  thermal  bounday  layer  region 


next  to  the  ceiling.  Incidentally,  the  isotherm  and  velocity  vector  scales 
in  Figs.  8  to  12  are  not  uniform,  and  consequently  the  diagrams  should  not  be 
used  quantitatively  for  comparison  purposes. 

CONCLUDING  REMARKS 

The  present  study  develops  a  three-dimensional  differential  field  model  based 
on  general  orthogonal  coordinate  systems  which  includes  the  physical  effects  of 
turbulence,  strong  buoyancy,  full  compressibility,  surface-surface  and  surface- 
flame  radiation  exchange,  pressure  rise  due  to  fire  loading,  and  heat  losses 
through  the  wall.  The  model  is  capable  of  dealing  with  geoometries  that 
represent  enclosures  with  combinations  of  shapes  which  are  describable  by 
orthogonal  coordinates.  The  model,  which  is  also  based  on  a  control -volume 
finite-difference  approach  with  primitive  variables,  staggered  cells,  and  QUICK 
scheme,  and  both  local  and  global  pressure  corrections,  is  used  to  simulate  a 
set  of  full-scale  fire  test  data  obtained  in  the  NRL  FIRE  I  fire  test  facility, 
which  is  a  pressure  vessel  with  a  cylindrical  mid-section  and  hemispherical  end 
caps  on  both  ends.  Despite  some  uncertainty  in  the  heat  release-rate  input 
data,  reasonable  comparisons  in  both  the  temperature  levels  at  selected  points 
in  the  test  tank  and  the  tank  pressure  between  the  calculated  and  test  results 
have  been  obtained.  However,  there  is  also  indication  that  the  calculated  tem¬ 
peratures  close  to  the  tank  ceiling  overpredict  the  test  values,  and  possible 
reasons  for  this  discrepancy  have  been  noted. 

Since  the  completion  of  this  study,  two  additional  capabilities  of  the 
field  model  have  been  developed  and  these  Include  effects  of  recirculating 
fans  to  provide  ventilation  inside  the  test  tank  and  the  placement  of 
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horizontal  deck  Inside  the  tank.  At  the  present  time,  the  field  model  treats 
the  fire  as  a  volumetric  heat  source  with  prescribed  heat  release  rate. 
Efforts  have  been  initiated  to  develop  a  combustion  model  so  that  the  fire 
envelope  and  heat  release  rate  do  not  need  to  be  prescribed.  Furthermore, 
plans  are  also  being  made  to  incorporate  a  gas  and  soot  (smoke)  radiation 
model  and  a  more  refined  turbulence  model  into  the  overall  field  model. 
Results  of  these  efforts  will  be  reported  at  future  dates. 
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mean  isobaric  heat  capacity,  Eq.  (4),  kJ/kgK 

radiation  view  factor 


components  of  gravity  vector,  m/s2 

radiation  matrix,  Eq.  (41) 

determinant  of  g..j 

covariant  metric  tensor 

contravariant  metric  tensor 

gravity  vector 

base  vector 

enclosure  height,  m 

coefficient  of  heat  transfer,  w/m2K 

scale  factor 

total  heat  flux,  w/m2 

constant  in  turbulence  model 

thermal  conductivity,  W/mK 

thickness  of  wall,  m 

mixing  length,  Eq.  (21),  m 

total  momentum  flux,  N/m2 

momentum  flux,  Eq.  (29),  N/m2 
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turbulent  Prandtl  number 

static  pressure,  N/m2 

heat  source  w/m3 

conduction  flux,  w/m2 

radiation  flux,  w/m2 

gas  constant,  Nm/kgK 

gradient  Richardson  number,  Eq.  (20) 

line  between  surfaces  i  and  j 

radius  variable,  m 

source  term,  N 

source  term,  Eq.  (31),N 

absolute  temperature,  K 

outside  ambient  temperature,  K 

velocity  components  in  system,  m/s 

velocity  components  in  01  system,  m/s 

matrix  components,  Eq.  (37) 

i nverse  of  X  .  . 

Cartesian  coordinates,  m 

angle  between  normal  to  A.  and  R..,  rad 

i  i  j 

control  volume,  m3 
Kronecke  delta 
surface  emissivity 
circumferential  angle,  rad 
orthogonal  coordinate 


viscosity,  kg/ms 


a 


Stefan-Boltzmann  constant,  w/m2K4 
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tensor , 
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1 J 

shear 
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in  0  system,  N/m3 

a  J 
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tensor , 

Eq.  ,28),  N/m2 

polar 

angle, 

rad 

Subscripts 

eff  effective  (laminar  +  turbulent) 

R  reference  quantities 

s  solid 

t  time  derivative 
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